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Abstract 

Recent Monte Carlo simulations of the critical point of the restricted primitive model for ionic 
solutions are reported. Only the continuum version of the model is considered. A finite size scaling 
analysis based in the Bruce- Wilding procedure gives critical exponents in agreement with those of 
the three-dimensional Ising universality class. An anomaly in the scaling of the specific heat with 
system size is pointed out. 

PACS numbers: 
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I. INTRODUCTION 



The primitive model (PM) for electrolytes, molten salts, colloids, etc is a mixture of M 
species of charged hard spheres living either on a lattice or within a continuous volume of 
real space. In this paper we shall focus only on the off-lattice version of the model. The 
simplest version of the PM consists of a binary mixture (i.e. M = 2) of positive and negative 
charged hard spheres ±q all with the same diameter a. Under this form the model which is 
thought to be the prototype of many ionic fluids has been christened the restricted primitive 
model (RPM). A thermodynamic state of the RPM is entirely specified by a reduced density 
p* = Na 3 /V (N number of ions, V volume) and a reduced temperature T* = kTa/q 2 (k 
Boltzmann's constant). 

The RPM undergoes a liquid-vapor transition which has been studied extensively these 
last past years by means of Monte Carlo (MC) simulations and various theoretical ap- 
proaches. In particular the behavior of the system at its critical point (CP) has been the 
subject of a huge amount of numerical and theoretical studies. The question is obviously of 
great importance since it is reasonable to assume that real electrolytes - or at last a large 
class of them- and the RPM belong to the same universality class which dictates a similar 
critical behavior. 

It is perhaps the right place to note that the main feature of ionic solutions is that the 
pair potential between two ions i and j at a distance = |r^ — r 3 -| which reads as 



v(rij) = — for r i:j > o 

= +00 for Tij < a (1.1) 

is a long range interaction. This fact would suggest classical (mean field) behavior, whereas 
the well-known screening of the interactions pleads in favor of an Ising-like criticality typical 
of systems with short range interactions. 

On the experimental side it seems well established now that for many real electrolytes 
apparent mean field behavior applies with sharp crossover (much sharper than in nonionic 
fluids) to Ising criticality close to the critical temperature^. 

At the moment there exists no convincing theoretical proof showing that the RPM be- 
longs to the Ising universality class^^ and only sophisticated Monte Carlo simulations can 



support this claim. Most numerical studies of the CP of the off-lattice version of the RPM 
were performed by the Orsay group and I would like to review our contributions towards a 
better understanding of the critical properties of this model in the lines below. 

II. A BRIEF HISTORICAL SURVEY 

Quite generally, a single component fluid will undergo a liquid vapor transition if the pair 
potential which is assumed to represent the molecular interactions is (sufficiently) attractive 
at large distances. From this point of view the situation is not so clear in the case of the 
RPM (cf eq. fjl.ljO and the very existence of the transition is not guaranteed. Several studies 
were necessary to clarify this point and a brief historical survey is worthwhile. 

The first evidence that the RPM actually undergoes a liquid-vapor transition can be 
tracked back to two papers of Chasovkikh and Vorontsov-Vel' Yaminov (CVVY) published 
as soon as in 1976 6,7 . These authors performed isobaric MC simulations and found a tran- 
sition with a CP located at T* = 0.095, p* = 0.24. Several years after (in 1991) Valleau 
studied three isotherms of the RPM with his method of the density scaling MC and ob- 
tained a different location for the CP, namely T* = 0.07, p* = 0.07 s . Subsequently (in 1992) 
Panagiotopoulos 9 obtained still different results, i.e. T* = 0.056, p* = 0.04, by performing 
MC simulations in the Gibbs ensemble (GE), at the moment a powerful new method of 
simulation which he had invented a little bit earlier—. Subsequent GE simulations using an 
improved biased MC samplingii yielded Panagiotopoulos and Orkoulas to the new estimate 
T* = 0.053, p* = 0.025. Finally, making use myself of the Gibbs ensemble combined with 
the use of hyperspherical geometries I obtained rather T* = 0.057, p* = 0.0 l 12 i 13 . 

Commenting on this striking dispersion of the MC data Prof M. Fisher talked once of the 
"sad street of numerical simulations". This was in 1999, at the SCCS conference, St Malo, 
France and, at this point of the story, I must agree with him retrospectively. However many 
advances have been done since. Before giving an account of these new achievements some 
comments are in order. 

• (i) All the MC studies confirm the existence of a liquid vapor transition for the RPM. 
It seems to take place at unusually low densities and temperatures. Caillol and Weis 
give further support for such a low critical temperature^. Moreover it turns out that 
the coexistence curve is very dissymmetric ^ 1 ^ 12 . 



• (ii) The MC simulation of ionic systems is a numerical challenge due to the long range 
of Coulomb potential. In order to deal with this, some caution is needed. Thus, in the 
case of MC simulations performed in a cubic box with periodic boundary conditions 
(PBC), one must use Ewald potentials in order to obtain the correct physic o 15 ' 16 ' 17 . 
The point is that the Ewald potential is the solution of Poisson equation in a cubico- 
periodical geometry^ 7 - and many properties of ionic fluids (electro-neutrality, screening, 
etc) are a consequence of this fact. In their MC simulations CVVY and Valleau 
considered truncated Coulomb potentials and very small samples of N = 32 particles 
which yields quantitatively wrong results. By contrast the data of Panagiotopoulos 
et a/j2»ii are more reliable since Ewald sums have been used. The same remark apply 
to my simulations which were performed on a 4D sphere (a hypersphere for short) by 
considering interactions obtained by solving Poisson equation in this geometry. This 
alternative method of simulation is therefore also indisputably correct, moreover it 
is much more efficient. The rough agreement observed between the simulations of 
refs.-^tll and^ 2 - both involving the same number of ions, i.e. N = 512, is therefore not 
fortuitous. 

• (iii) None of the above mentioned studies took correctly into account finite size ef- 
fects which are of an overwhelming importance near a CP. These effects affect the 
behavior of finite systems as soon as the correlation length of the critical density fluc- 
tuations is of the same order of magnitude than the size of the simulation box. In the 
simulation o 9 i n i 12 some "apparent" critical temperature T* has been measured which 
could be very different from its infinite volume limit T c *(oo). 

In order to extract from MC simulations the critical behavior of the RPM in the thermo- 
dynamic limit (i.e. the critical exponents) and also the infinite volume limit of T* and p* it 
is necessary to perform an analysis of the MC data in the framework of the finite size scaling 
(fss) theory which is part of the renormalization group (RG) theor y 1 ?' 19 . In this approach 
one needs to work in the Grand Canonical (GC) ensemble rather than in the Gibbs ensemble 
which is ill adapted for a fss analysis. Subsequent MC simulations on the RPM were thus 
all performed in this ensemble. Panagiotopoulos and coworkers turned their attention to 
the lattice version of the RPM whereas the Orsay group continued to work on its off-lattice 
version. 



III. FINITE SIZE SCALING ANALYSIS OF MC DATA 



A. Scaling fields and operators 

Starting with the seminal work of Bruce and Wilding (BW) 20 i 21 i 22 simulation results for 
the critical behavior of fluids have customarily been analyzed along the lines of the so-called 
revised scaling theory of Rehr and Mermin 2 ^. In this approach one first defines scaling fields 
and operators aimed at restoring the particle-hole symmetry and therefore to map the the 
fluid onto a magnetic sytem with Ising-like symmetry. 

The two relevant scaling fields h (the strong ordering field) and r (the weak thermal 
field) are assumed to be linear combinations of deviations from their critical values of the 
chemical potential /i and the inverse temperature (3 = 1/T (reduced values are assumed 
henceforward). One thus has 



h = /I - ii c + r{j3 - Pc) 

t = (3 C - + s(fj, - fji c ) , (3.1) 

where r and s are the field mixing parameters which define the mapping. Of course relations 
()3.1|) are valid only in the vicinity of the CP. The conjugate scaling operators Ai and £ are 
then defined as 

1 d 1 
< M > = TTTTT-lnS = (< p> -s <u>) 

V oh 1 — sr 

X c) 1 

<£> = —— lnS=- {<u>-r<p>), (3.2) 

V or 1 — sr 

where S is the GC partition function of the RPM, p the total number density, and u the 
internal energy per unit volume. Brackets < . . . > denote GC averages. M. is the order 
parameter (magnetization) of the magnetic system associated with the fluid and £ its mag- 
netic energy. £ should be invariant under the transformations (Ai,h) — > (—Ai,—h) for 
appropriate values of s and r. In this framework the coexistence curve is therefore defined 
by the eq. h = 0. 

The revised scaling of Rehr and Mermin implies the analyticity of the coexistence chemical 
potential /x(T) at T*. Although this is the case for some peculiar lattice gas models with 
"hidden" symmetries there is no reason that in general, for fluid systems fi(T) should lack a 



singularity as recognized already by Rehr and Mermin^ 3 . and emphasized more recently by 
Fisher and co-worker a 40 ' 41 ' 42 . 

B. The scaling hypothesis 

A central role in the subsequent fss analysis is played by the GC joint distribution 
Vl{M,£) oc Vl{p,u) for the scaling operators M. and 8. Following B W 20 ' 21 ' 22 we will 
assume that, in the immediate vicinity of the CP, Vl(M,£) obeys to the following scaling 
law: 

V L (M,8) = a^af L d - y - L d - yh V(a^L d - yh 5M, 

. . . a^L d - yT 58, a M L yh h, cl £ L Vt r, ai L Vi ,...), (3.3) 

where L are the linear dimensions of the system (taken as V 1 ^ 3 , where V is the volume of 
the simulation box, either a cube or a hypersphere). I have denoted by SAi = M.— < M. > c 
and 58 = 8— < 8 > c the deviations of the scaling operators from their value at criticality. 
The cornerstone of this scaling hypothesis is that the function V which enters eq. (J3.3J) is 
universal in the sense that it depends only upon the universality class of the model and of the 
type of geometry considered. The constants a_M, a>£, and aj are system dependent constants 
which are defined in such a way that V has unit variance. Finally, the renormalization 
exponents i/h, y T , and y,i which enter eq. are defined as 

y h = d — P/v 
Vt = 

Vi = -9/v (3.4) 
in terms of the usual critical exponents: 

• (3 exponent of the ordering field, i.e. < 5A4 >~ \t\@ for T* < T* at h = 

• v exponent of the correlation length, i.e. £ ~ |r| ~ v 

• 9 Wegner's correction-to-scaling exponent (first irrelevant exponent). 



The scaling hypothesis (|3.3|) was established on a solid RG basis for Ising-like systems^ and 
received substantial supports from MC studies — . We stress once again that the coexistence 
curve is determined in this approach by the condition h = and that, at coexistence, the 
order parameter distribution Vl{M) should be an even function of M.. In practice this 
symmetry requirement can be satisfied by tuning the two parameters (/i, s) at a given {3. 
We now concentrate our attention on the scaling behavior of the histogram Vl{M)- 



C. The matching procedure 

Integrating both sides of eq. ()3.3j) over £ one finds that, along the coexistence line h = 
one has 

V L (M) = ajj L d - yh V (a^L d ~ yh 5M, a e L* r, a^) , (3.5) 

where, in the r.h.s. the dependence of the universal function V upon h has been discarded 
for clarity Let us define now x = aJ^L d ~ Vh SAi, then, assuming r ~ and L ~ oo a Taylor 
expansion of eq. ()3.5|) yields 

V L (M) = a]^ L d - yh \v*(x) + a £ D>T r V\{x) 

+ afL 2yT t 2 V* 2 (x) + ... + ai L yi Vl{x) + ...1 , (3.6) 

where the various V* entering the r.h.s. are universal functions. Note that, for L = oo the 
normalized ordering field distribution Vl{M) collapses onto an universal function V*(x) at 
t = 0. For L finite but large Vl(M) collapses approximately onto V*(x) at some apparent 
tl oc L~ VT+Vi . Since for h = one has r oc f3 — (3 c then the matching of the histogram Vl{M) 
onto the universal function V*(x) should occur at some apparent temperature T*(L) scaling 
with system size as 

T*(oo) - T*{L) oc L- {e+1)/u + .... , (3.7) 
where T*(oo) denotes the infinite volume limit of the critical temperature. 



D. Technical details 



To assess the critical behavior and the critical parameters of the system, we need, in a 
first step, to locate the coexistence curve h = 0. At a given temperature j3 close to (3 C the 



ordering distribution function Vl{M) depends solely on the chemical potential \x and the 
mixing parameter s. At coexistence, the value of (fi,s) can be obtained unambiguously by 
requiring that Vl{M) is symmetric in M.— < M. > 22 . Tuning at will \x and s at given f3 
requires to know the joint histogram Vl(M,£) oc Vl{p, u) for a continuous set of values of 
fi at a given {3. Moreover, since this analysis must be performed at different (3 one needs in 
fact to know Vl{M,£) for a continuous set of values of in the critical region. This 

technical difficulty is circumvented by using the multiple histogram reweighting proposed by 
Ferrenberg and Swense n 26 i 27 i 28 . With this method one can obtains ' M (p, u) for a continuous 
set of values of from the knowledge of R histograms Pf" M '(p, u), i = 1, . . . , R obtained 
by performing R distinct MC simulations in the R (neighbor) thermodynamic states (/3j, /ij). 

Since the precision of the simulations of fluid systems has still not reached that obtained 
in the MC simulations of Ising like systems it is impossible to construct ex nihilo the fixed 
point universal distribution V*{x). In refs . 29 ^ 30 our attempts to match Vl{M) on V*(x) 
were realized by using the estimate of V* s (x) made by Hilfer and Wilding 32 for the 3D Ising 
model. Two new -and better- estimates of V* s (x) obtained by Tsypin and Blote 33 - for the 
3D Ising model and the spin-1 Blume-Capel model were considered in ref.— . The discussion 
is postponed to next section. 

IV. RESULTS 

A. General discussion 

It turns out that the field mixing parameter s of the RPM is practically independent of 
the temperature and of the size L of the system. Its magnitude, s ~ — 1.46 29 i 30 i 31 , is much 
higher than for neutral fluids (typically s ~ 0.02 for square well or Lennard- Jones fluids*^) 
which explains the large dissymmetry of the liquid-vapor coexistence curve of the RPM. 

The collapse of the ordering operator distribution Vl(-M) onto the universal ordering 
distribution V* s (x) given by the Blume-Capel model^ 3 - is depicted in Fig. 1 for four different 
values of the volume ranging from V/a 3 = 5000 to V/a 3 = 40000, i.e. up to a linear size 
L/a = 34. At volume V/a 3 = 5000 a mismatch is observed at the lowest values of Ai due 
to an inadequate sampling of of the low density configurations at small volume. The overall 
good agreement leads us to conclude that the universality class of the RPM is that of the 



3D Ising model. 

The reduced apparent critical temperature T*(L) versus the size L of the system (in 
reduced units) has been plotted in Fig. 2. Depending on the choice made for the universal 
ordering distribution V* s {x) one obtains two sets of values of T*(L) from which T c *(oo) can 
be obtained by making use of eq. (|3.7j) . One obtains T c *(oo) = 0.04917 ± 0.00002 using 
V* s {x) derived from the Blume-Capel model and T c *(oo) = 0.04916 ± 0.00002 using V* s (x) 
obtained for the 3D-Ising model. The approximate V* s (x) of Hilfer and Wilding yields 
slightly different results. Note that in all cases we have used the Ising values v = 0.63C 35 
and 9 = 0.53 36 of the critical exponents. 

The previous analysis merely establishes the compatibility of the MC data with an Ising- 
like criticality. One can try to go beyond by considering the scaling behavior of the Binder 
cumulant 

As a consequence of the scaling hypothesis ()3.6|) one can show that, at coexistence (h — 0), 
Qb{L) should scale with system size as 



Qb{L) = Q c + q l ((3-(3 c )L 1 ^ + q 2 ((3-f3 c ) 2 L 2 / v 

+ q3(P-Pc) 3 L 3 ^... + hL yi + ... (4.2) 

where the last term takes into account contributions from irrelevant fields and qi, q 2 , ?3, 
and bi are non universal constants. If the contribution of irrelevant fields could be neglected 
then the curves Qb{L) would intersect at the fixed point Q c . As apparent in Fig. 3 this is 
clearly not the case and corrections to scaling must be taken into account. 

Recall that for the 3D-Ising model the fixed point value is Q c = 0.623^ and that the 
exponent of the correlation length has the value v = 0.63C 35 . We have attempted to fit all 
our MC data with eq. (j4.2j) . If all the parameters in the RHS of eq. ()4.2|) are kept free such 
an ambitious fit turns out to be impossible. Various other less satisfactory strategies can be 
considered however. 

• Fixing (3 C = 1/0.04917, yi = —Qjv = —0.84 and leaving free all the other parameters 
one finds a fit better than 1 per cent and Q c = 0.63 ± 0.01, v ~ 0.66 ± 0.03. 



• Conversely, fixing Q c = 0.623 and 9 = 0.53 one obtains (3 C = 1/0.04918 and v ~ 
0.63 ±0.03. 

The variations of Qb{L) as a function of (3 for the different volumes is shown in Fig. 3. 
Although there is considerable spread in the intersection points due to correction-to-scaling 
contributions, the corresponding values of Q c are close to the Ising value Q c = 0.623 and 
permit to rule out mean field behavior (i.e. Q c = 0.457—). 

Further support for Ising criticality is provided by the behavior of < 5M 2 > at T*(L). 
According to the scaling hypothesis (J3.6|) it should scales as L 2/3 ^ u with system size. From 
the slope of the curve displayed in Fig. 4 one obtains (3/v = 0.52 in accord with the 3D Ising 
value (0.517) and in clear contrast with the classical value 1. 

In summary, our fss analysis leads to an estimate of the critical exponents v and f3/v and 
the Binder cumulant Q c based on the sole knowledge of the critical temperature and the 
renormalization exponent 9. Within the numerical uncertainties these values are compatible 
with Ising-like criticality. Our conclusion is that the RPM, as ordinary neutral fluids, belongs 
to the universality class of the Ising model. 

A complete discussion of our MC data is out of the scope of the present paper and 
can be found in ref.— . For completeness I give below the values obtained for the critical 
temperature, chemical potentials and densities (the critical pressure is largely unknown): 

• T* = 0.04917 ±0.00002 

• p* = 0.080 ± 0.005 

• li* = -13.600 ±0.005 

B. The specific heat 

The revised scaling theory of Rehr and Mermin which is the framework of our fss analysis 
is however not the most general scaling theory which can be proposed for a fluid system 
lacking the "particle-hole" symmetry. Its main weakness, as recognized already by Rehr and 
Mermin^, Yang and Yang22, and more recently by Fisher and co-workers^*^^, is that it 
assumes the analyticity of the chemical potential at coexistence p(T) at the critical point. 
The more general scaling assumption should yield singularities for both p(T) and p(T) as 



T* — > T*. Let us examine the consequences of these singularities on the behavior of the 
specific heat capacity at constant volume Cy. In the two phase region it can be rewritten 
as^ 



Cy 




V 



-NT 



or 2 




V 



(4.3) 



where C p (not to be confused with the specific heat capacity at constant pressure) and C M 
(not to be confused with the specific heat capacity at constant chemical potential) denote 
the two contributions to Cy. I stress that, in eq. ()4.3|) p(T) and p(T) denote the pressure 
and the chemical potential at coexistence. The formula can be used for any density p g {T) < 
p < pi(T) within the two phase region (p g (T) and pi(T) being the densities of the gas and 
the liquid at coexistence respectively). In the revised scaling theory only C p diverges as 
\T* — T*\~ a whereas one expects a divergence of both C p and (both as |T* — T*\~ a ). In 
Fig. 5 I display the curves C M (T) and Cy{T) along the locus 



for the four volumes considered in our last MC simulations 31 . Although the peak posi- 
tions shift correctly as oc L~ x l v with system size, in accord with fss theor y 1 ?' 1 ? , there is no 
detectable scaling of the heights of the peaks which should scale as L a l v with L. These ob- 
servations corroborate similar results obtained by Valleau and Torri o 43 ' 44 . In particular C M 
does not show any anomaly which should challenge the use of eqs. (|3.1|) for the scaling fields. 
A possible explanation for the non singular behavior of Cy(T) is that the amplitude of the 
singular term in Cy(T) is small in the RPM and that the specific heat is dominated by its 
regular part. Note however that, the peak heights in Cy(T)/V would scale, assuming Ising 
value for a only by a factor 2 a ^ u ~ 1.12 when doubling the linear dimensions of the system. 
It is possible that such a small effect is not detectable within the statistical uncertainties of 
our calculations. 

V. CONCLUSION 



Xnnn oc< (N- < N >) 3 >= , 



(4.4) 



In this paper which resumes my talk at the Lviv NATO workshop I have described 
recent attempts to elucidate the nature of the critical behavior of the RPM model for ionic 



fluids, prototype of a system governed by long range Coulomb interactions by means of MC 
simulations. After endeavor over more than a decade we have now reached a point where 
we can claim confidently that the RPM belongs to the same universality class than the 3D 
Ising model. The critical values of non-universal quantities such as the temperature and the 
chemical potential were established with a high accuracy whereas the uncertainties on the 
critical density are more significant, and the critical pressure is unknown. 

The behavior of the constant volume specific heat gives no indication of the expected 
L"/" scaling of the peak height within the range of system sizes considered in the most 
recent simulations. Recent investigations of Camp and co-workers 45 where differences in the 
behavior of Cy in the canonical and the GC ensemble are reported have emphasized this 
problem. At the moment it is difficult to explain this unexpected behavior of the specific 
heat. 

I have only discussed the properties of the continuous version of the RPM. The phase 
diagram of the various lattice versions of the model is in fact more comple x 46 ! 47 and was not 
described here due to a lack of place. I have also excluded from my presentation assymetric 
versions, either in charge or/and in size, of the continuum or lattice versions of the primitive 
model. The interested reader should consult recent works of Panagiotopoulos et q/ , 48 i 49 i 50 
and de Pablo et a / . 51 i 52 i 53 . 
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ure Captions 

• figure 1: Collapse of the ordering distribution Vl{M) onto the universal Ising order- 
ing distribution V* s (x) for V/a 3 = 5000, T*(L) = 0.004934, s = -1.45; V/a 3 = 10000, 
T*(L) = 0.004926, s = -1.465; V/a 3 = 20000, T*(L) = 0.004921, s = -1.47; and 
V/a 3 = 40000, T*{L) = 0.004922, s = -1.43. V* s (x) (solid circles) if the MC result of 
Tsypin and Blote (Ref.—) obtained for the Blume-Capel model. The scaling variable 
is x = aJ^L^^S Ai where is chosen in such a way that has unit variance. 

• figure 2: The apparent critical temperature T*(L) function of L-^ e+1 ^ u with 9 = 
0.53, v = 0.630 obtained by matching the universal ordering distribution calculated 
for the Blume-Capel model (top) and the Ising model (bottom). Extrapolating by 
linear least square fit to the infinite volume limit yields T*(oo) = 0.04917 ± 0.00002 
(top) and T c *(oo) = 0.04916 ± 0.00002 (bottom). 

• figure 3: Variation of Qb{L) as a function of the inverse temperature (3 for the 
different volumes considered in ref.—. From top to bottom V/a 3 = 40000, 20000, 10000, 
and 5000 respectively. The symbols are the MC data and the lines the fits obtained 
by means of eq. ()4.2|) . 

• figure 4: Variations of In < 5A4 2 > at T*(L) as a function of InL. The slope of the 
linear least square fit is 2(3 /v = 1.04. 

• figure 5: Variations of the total specific heat at constant volume Cy/V and the 
contribution C^/V with temperature along the locus Xnnn = at volumes V/a 3 = 
5000, 10000, 20000, and 40000 (from left to right). 
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